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1. Introduction 



The square double quantum well (DQW) often is used as a toy model to demonstrate 
the interaction between quantized energy levels due to particle tunneling through a 
potential barrier separating individual wells [U [2} EJ HJ |5j [6] . Recently the DQW model 
had attracted considerable attention in semiconductor heterostructure physics because 
of its applications in nanoelectronics [7J [8j |9]. The tunneling conductance properties 
of semiconducting DQW devices as well as drag effects that result from interaction 
between electrons moving at different velocities in different wells was recently discussed, 
for example, in review articles [TOj [TTj. 

Appearance of transcendental equations that describe DQW spectrum limits direct 
application of analytical methods in tackling the eigenfunction problems. Initially the 
problem of finding the eigenfunctions has been solved by perturbation theory assuming 
that energy level splitting due to tunneling is small The most recent analytical 
approach heavily relies on symmetry properties of the DQW [6j. Of course, this 
restriction can be relaxed by resorting to numerical methods J2j HJ [6j EJ [12]. However, in 
many cases a knowledge of analytical form of the wave function is more preferable. For 
example, in the wave packet dynamics problems the closed form solution allows one to 
construct a direct superposition of eigenfunctions to make a computational task easy. 
Here we demonstrate that one can push the problem further and calculate the relevant 
eigenfunctions exactly by exploiting a computer based Grobner basis algorithm jT3j . In 
sections [2] and |3] the spectrum and eigenfunctions of a general DQW are calculated using 
the Grobner basis, and in section H] the results are applied to find closed form expression 
for optical dipole matrix element of the DQW. 



2. Spectrum 



The one-dimensional DQW with flat potentials in each of regions 1 — 5, as shown in 
figure [H is described by the following piecewise function of coordinate x 



V(x) 



V c if x < 

if < x < a 

Vb ii a < x < a + b 

if a + b<x <2a + b 

V c if x > 2a + b 



(2.1) 



where V c is the confining potential (referenced from the bottom of wells) and V& is the 
height of central barrier separating two identical quantum wells. The mirror symmetry 
of the system ensures that the quantum states of such a DQW have either even or odd 
parity. 
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Figure 1. Symmetric double quantum well with central barrier of width b and height 
Vf). The eigenenergy E n is referenced from the bottom of wells of width a. Electron 
effective mass in regions 1 — 5 is assumed to be different. 



Only bound states will be considered here. These states can be normalized to unity 
over entire x axis. The wave function ifj(x) in the regions 1 — 5 has the following shapes: 

fa —A\ sin kx + C\ cos kx, 

fa =B 2 e Xb(a - x) + B 3 e- Xb{a+b - x) , (2.2) 
^4 =^2 sin k(2a + b — x) + C-i cos k(2a + b — x) , 
fa =B i e Xc(2a+b - x) , 



where k is the free-electron wave vector, k = ^ 2m$E / h 2 , in the quantum wells 
of width a. The energy E is referenced from the bottom of the wells. The wave 
vectors of evanescent waves in the exponents are Xb = a/ (2m b /h 2 )(Vi > — E) and 
Xc = a/ (2m c /h 2 )(V c — E), where we have introduced different electron masses, namely, 
mo inside the wells, mj in the barrier and m c in the confining potential. This is typical 
to semiconductor heterostructures, where the DQW is made of nanometer layers having 
different forbidden energy gaps. As a result, the electron effective mass depends on 
coordinate x. 

In equations (12. 2p there are eight unknown coefficients that must be calculated. 
Because of symmetry, the number of coefficients, in fact, can be reduced. However 
we shall not do this since the Grobner basis algorithm will take account of symmetry 
properties of polynomials automatically. The standard BenDaniel-Duke boundary 
condition [14] which takes into account mass difference on right (r) and left (I) sides of 
the potential step at coordinates X = 0, X = a, X = a + b and X = 2a + b will be used 

A(X + ) =ih{X~), (2.3a) 



1 dlp r 

m r dx 



x+ mi dx 



(2.3b) 
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(2.6) 
(2.7) 



Equations ( 12. 2 p and the boundary conditions yield the system of eight linearly dependent 
equations 

B 1 -C x = 0, -A^/mo + B lXc /m c = 0, (2.4a) 

-B 2 - B 3 e~ bXb + A x sin ak + d cos ak = 0, (2.46) 

(B 2 Xb — B 3 Xbe~ bXb ) /m b + (Aik cos ak — Ciksmak)/m = 0, (2.4c) 

B 3 + B 2 e~ bXb - A 2 sin ak - C 2 cos ak = 0, (2 Ad) 

(B 3 Xb — B 2 x b e~ bXb ) / m b + (A 2 k cos ak — C 2 k sin ak)/m = 0, (2.4 e) 

-B A + C 2 = 0, -A 2 k/m + B 4X c/m c = 0. (2.4/) 

The determinant D that follows from this system determines the spectrum of discrete 
energy levels of DQW. The symmetry of the problem ensures the factorization of the 
determinant 

D = -mQ 4 m- 2 m^ 2 e- 2bXb D s D a = , (2.5) 
where D s and D a refer, respectively, to symmetric and antisymmetric states, 
D s = - km [(xcm b - x 6 m c ) + e bXb (x c m b + Xbm c )] cos ak+ 
\{k 2 m b m c + XbXcm 2 ) + e bXb (k 2 m b m c - XbX c ml)} sin ak , 

D a =km [(x c m b - Xbm c ) - e bXb (x c m b + Xbm c )] cos ak- 

[(k 2 m b m c + XbXcml) - e bXb (k 2 m b m c - XbXcml)] sin ak . 

To advance further the transcendental equations D s (k) = and D a (k) = which 
determine, in turn, the spectrum of symmetric and antisymmetric discrete energy levels 
have to be solved explicitly. Unfortunately these transcendental equation only can be 
solved by numerical methods. If DQW parameter values are known, then roots of (12.61) 
and (12. 7p define the spectrum of all wave vectors k n , or equivalently discrete eigenenergies 
E n = h 2 k 2 l /2m of the DQW, where n is the energy level index. 

In a special case when the DQW heterostructure is fabricated from two types of 
nanolayers (labelled b and 0) we have that V c = V b and m c = m b . Then Xc — Xb, and 
the determinants ( 12. 61) and ( 12. 71) simplify to 

D Sja = — 2A;x^ ) mom^ ) e bXi, cos ak± 

[(k 2 m 2 b + xWo) + e bXb (k 2 m 2 b - xW,)} sin ak = , 2 ' § 

where plus/minus signs correspond to symmetric/antisymmetric states. When m = m b 
further simplification is possible 

D S:a = 2 cos ak + (f - sin(a£;) ± (f + sin(aA;)e- x6 = 0, (2.9) 

where now k = ^2^E/h, x = ^2m (V - E)/h and f = x/k = ^{V -E)/E. 
Here the plus/minus sign corresponds to the antisymmetric/symmetric state relative to 
the center of the DQW structure, respectively. The expression (12.91) can be found in 
references [3, [12] , where the energy in the presented formulae is counted from the top of 
the wells. When the barrier width b — > oo, equation (12. 9p goes back to the well known 
formula for an isolated quantum well. 
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When the particle energy E is larger than the height V b of the barrier but smaller 
than the confining potential, Vj, < E < V c , the particle still remains localized. The only 
difference is that in the regions 2 — 4 wave function now oscillates, i.e. the eigenf unctions 
ip{x) here are described by trigonometric functions only. It is easy to see that the above 
solution at E < V b remains valid if we account for hyperbolic functions properties 
sinh(ix 2 ) = isinx 2 , cosh(ix 2 ) = cos X2 and notice that in this case Xb can be replaced 
by i X b = i^(2m /h 2 )(E-V b ). 

3. Eigenfunctions 

The coefficients in the wave function (12.21) depend on k n . Since the spectrum k n (or 
E n = h 2 k 2 /2m Q ) is determined by roots of the transcendental equations (12. 6p and 
(12. 7p . one is obliged to solve these equations using numerical methods. Nonetheless, 
as we shall see, the eigenfunctions can be explicitly calculated with the help of Grobner 
basis algorithm [13j [15] without any reference to the roots at all. Roughly speaking, 
a Grobner basis for a system of polynomial equations is a different system of simpler 
polynomials having the same roots as the original ones. Calculation of the Grobner 
basis to some extent resembles reduction of square matrix to triangular matrix. For 
further calculations it is convenient to introduce the following half angle variables 

x = tau(Wfe/2), y = tau(afc/2) (3.1) 



and express sine and cosine functions in (\2Aa\i-(2Aj) and (12. 6 p (or (12. 7p in case of 



antisymmetric eigenfunctions) through polynomial variables x and y, 

, 2x , 1 — x 2 

sin a/c = -, cosafc 



1 + X 2 1 + X 

sin bk = — —, - , cos bk J 



(3.2) 



1 + y 2 ' 1 + y 2 ' 

Calculating Grobner basis for coefficients A, B and C and requesting that new variables 
x and y to be eliminated, the Mathematica program generates basis which consists of 
146 polynomials. However, it should be stressed that the program can find the Grobner 
basis only if the spectrum equation, either (12. 6 j) or ( 12. 7ft is appended to the original 



polynomial system ( 12.4 a\i -( 2.4/). The following simplest polynomials were selected for 
symmetric states 



a - a - r Xc7Jl ° 

km c 

B\ = B± = C\ = C2S, 



(3.3) 



±C 2s m b e bXb ^Jk 2 m 2 + xl^l 



m c [k 2 m 2 (l + e b »)2 + x 2 b m 2 (-l + e b »)2]V2 ' 

where C 2 was replaced by C 2s to identify the state symmetry. The choice of 
sign of B2 and B% coefficients has to ensure derivative continuity at points a and 
a + b. It is straightforward to check that the solution ( 13. 3 j) indeed satisfies the 



initial equations (\2Aa\\-{2Aj). In ( 13. 3p all amplitudes are expressed through a single 
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coefficient C 2s} which in turn can be found from the normalization condition of the total 
wave function ip(x). The integration over x axis gives the normalization constant in the 
form 



C 2s = km c {G x + G 2 ) 



-1/2 



where 



Gi = Xc l [ k2m l( l + a Xc) + m x 2 c (m c + m ax c )} , 

m b (k 2 m 2 c + xl m l)[bXbk 2 m b + (k 2 m b + x b m o) smhbxb] 



G 



(3.4) 

(3.5) 
(3.6) 



Xb[k 2 m 2 b - xlml + {k 2 m\ + Xb m o) cosh b Xb] 

If all masses are assumed to be equal (mo = m& = m c = 1) the normalization constant 
simplifies to 



G 



2s 



k 2 



aXb + ^ + 

Xc 



Xl , b Xbk 2 + (k 2 + xl) sinh bxb 



k 2 ~ X\ + ( k2 + Xl) cosh&Xb' J 



-1/2 



.(3.7) 



Quite similar calculation for antisymmetric (C 2 — > C 2a ) states yields 

X c m 



B 2 



B 4 
B 3 



■C x = G, 



2a- 



-Ax = A 2 



C 



2a ■ 



km r . 



±C 2a m b e bXb {k 2 m 2 c + Xc m o) 1/2 



(3.8) 



m c [k 2 m 2 {-l + e fe »)2 + x 2 b m 2 {l + e b ^) 2 ] 1/2 ' 

where the choice of sign again follows from the derivative continuity condition. The 
normalization constant in this case is 

C 2a = km c {H l + H 2 )- 1 ' 2 , 



where 



H \ = X c X \k 2 m 2 c {\ + axe) + m x 2 c (m c + m ax c )] , 

m b (k 2 ml + x 2 c m l){ - bxbk 2 m b + (k 2 m b + Xb m o) sinh&Xfo] 



H- 



Xb[ ~ k 2 m 2 b + xlml + (k 2 m 2 b + xl m o) coshfexb] 
When all masses becomes equal the normalization constant C 2a reduces to 

Xl , -bXbk 2 + (k 2 + xl) sinh bx b 



(3.9) 

(3.10) 
(3.11) 



2a 



1 + 



k 2 



aXb^ 



X 2 C ~k 2 + xl + {k 2 + xl) cosh bxb 



-1/2 



.(3.12) 



As far as a more general non symmetric DQW problem concerns, the calculations of 
the Grobner basis indicates that, in contrast to solutions (I3.3P and (13. 8p . at least one 
of the coefficients A, B, or C includes the trigonometric functions. In this case the 
determinant D does not factorize to symmetric and asymmetric parts either. 
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b(nm) 



Figure 2. The energies of GaAs/Gao.gAlo^As DQW as a functions of the central 
barrier width at a = 6 nm. 



d r , 



4. Dipole matrix element 

The knowledge of eigenfunctions allows one to carry on with analytical calculations. As 
an example we shall find closed form expression for dipole matrix elements between even 
ip s (k n ,x) and odd ipa{k m ,x) discrete states 

ip*(k n ,x)lx - a- -jip a (k m ,x)dx = 2d! + 2d 2 + d 3 . (4.1) 

Here the subscripts s and a refer to, respectively, even and odd symmetry states and 
di is the contribution of the 2-th region indicated in the figure [TJ For a general case 
the expressions for dipole components d nStTna are rather complicated [16j. For simplicity 
below we present the expressions for the case when masses in all regions are equal, 
m c = rrib = mo and the central and confining barrier heights coincide, Xc = Xb = X- 
Since the energy of symmetric and antisymmetric states differ the wave vectors k and x 
are supplied by indices s and a. Thus the dipole expression have two kind of the wave 
vectors k s and k a , and evanescent modes Xs an d Xa- 

In the first and fifth regions the contribution to dipole is 

d\ = <i 5 = d\N / d\D , 



where 



dm 



k s k a [2 + (2a + b)(x s + Xa)}r s r a , 
VXs[(k 2 s - xj) + (k 2 s + xj) cosh bxs 
VXa[-(kl- xj) + Wl + xj) coshbxa] 



(4.2) 



(4.3) 



and 



diD 
d. 



6 a 



2(X, + Xa)W(k 2 + Xm 2 a + Xl)S s S a , 

- x 2 a (l + aXs) + k 2 s [l + (a + b) Xs ] + 



(4.4) 



(k 2 + X»)[(! + a Xs) cosh6x s + sinho^s 
X 2 a il + a Xa ) ~ k 2 a [l + (a + b) Xa \ + 



1/2 
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15 a) 



d(nm) 
10 
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6 
4 

2- 
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1 
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15 b) 

Figure 3. a) Dipole matrix elements <ii Si 2a and c?2a,3s as a function of barrier width b. 
b) Contribution of individual regions to dipole matrix d ls ^ a - 



♦ 


d\,d*, 




d2^4 




fi?3 




1/2 



)- 

In the second region it is 

d2 = d 2 N/d 2 D, (4.5) 



where 





= - ^r s r a [(k s - k a ) 2 (p 1 - 


Ps) - (k s - 


f k a ) 2 (p 2 -Pa)} , 


(4.6) 


Pi 


= [b(k s + k a )(k s Xa + kaXs) 


+ 2k s k a - 


2xsXa] cos a(k s + k a ) , 






= [b(k s - k a )(k s Xa ~ k a Xs) 


— 2k s k a — 


ZXsXa] cos a(k s - k a ) , 




P3 


= [b(k s + k a )(k s k a - XsXa) 


- 2k s Xa - 


2k a Xs] sin a(k s + fc a ) , 




P4 


= [b(k a - k s )(k s k a + XsXa) 


+ 2k a Xs - 


2k s Xa\ sin a(fc s - fc a ) , 




^2D 


_ {hi ~ Kf i 

{Xs + Xaf 






(4.7) 



One can see that trigonometric functions, which will give oscillations of matrix elements 
vs. the well width, appear only here. 

The third (barrier) region contribution to dipole is 

d 3 = d 3N /d 3D , (4.8) 
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where 



4e^* s+ * a )A; s A; a r s r> 1 cosh ^ + v 2 sinh , (4.9) 



2 2\ & Xs , _. 



and 



«i = - &Xa(x s - Xa) cosh — + A XsXa sinh o 
^2 = - 2(Xa + Xa) cosh ~Y + 46 X«(x5 ~ Xo) sinh ~7T > 

4d = ^(x* - Xa) 2 s s s a rf 1D , (4.10) 



;i + e^)2£;2 + (_l + e^)2 x 2 



-l + e fe ^) 2 A;2 + (1 + e b Xa y x 2 a 

k l + xl 

Figure [2] shows the dependencies of the first two energy levels E ls and E 2a as a 
function of the inner barrier width. The following parameter values that are typical 
to GaAs/Ga .8Al .2As DQW heterostructures, were used for production of pictures: 
a = 6 nm, b = (1 — 15) nm, V c = Vb = 0.1671 eV, m = 0.067m e , m c = m b = 0.0836m e , 
where m e is the electron mass in the vacuum. The increase of energy difference between 
levels with the decrease of b is assigned to tunnel coupling of levels. 

Figure [3^ demonstrates, respectively, the size of optical dipole matrix elements 
between a pairs of adjacent levels, c? ls 2 a and G?2 a ,3s> as a function of barrier width. 
Figure [3b shows the contribution of individual regions to the dipole di Sy 2 a - It is clear that 
a general trend and magnitude of dipole elements in figure [3^ can be understood if one 
assumes that only quantum wells contribute to the total dipole. In this approximation 
the functions ip\ = ip3 = V's — while the ip 2 and ip± can be approximated by half-period 
sine functions. Then du,2a reduces to 

, 2 f a . tix ( b\, . 7tx . a + b fA ,,\ 

d\s.2a ~ - / sm — \x-a-- (-sm — )dx = — — . (4.11) 
a J a V 2/ a 2 

The formula shows that dipole size increases linearly with the barrier width b as 
long as b remains much smaller than exciting light period. For 2a — 3s optical 
transitions one of sines should be replaced by sin(27rx/a). Then, similar calculation 
yields G?2a,3s ~ 16a /9ir 2 , which is independent of barrier width. The deviations from the 
obtained expressions in figure come from the evanescent mode contribution in barrier 
and confining potential regions. 

In conclusion, the presented example shows that application of Grobner basis 
algorithm in some cases allows to find closed form expressions for the total wave function 
and, therefore, to calculate the dipole matrix elements exactly without directly solving 
the transcendental equations that determines the spectrum of the DQW. Of course, 
the described method can be applied to other quantum systems for which eigenvalue 
equations cannot be explicitly solved as well. 
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